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Abstract 

In recent years there has been a growing interest in the fractional Fourier transform 
driven by its great number of applications. The literature in this field follows two main 
routes. On the one hand the applications fields where the ordinary Fourier transform can 
be applied are being revisited to use this intermediate time-frequency representation of 
signals; and on the other hand fast algorithms for numerical computation of the fractional 
Fourier transform are devised. In this paper we derive a Gaussian-like quadrature of 
the continuous fractional Fourier transform. This quadrature is given in terms of the 
Hermite polynomials and their zeros. By using some asymptotic formulae we are able to 
solve the quadrature by a diagonal congruence transformation equivalent to a chirp-FFT- 
chirp transformation, yielding a fast discretization of the fractional Fourier transform 
and its inverse in closed form. We extend the range of the fractional Fourier transform 
by considering arbitrary complex values inside the unitary circle and not only at the 
boundary. Interestingly enough, the congruence transformation evaluated at z — i, which 
gives the Fourier transform, improves the standard discrete Fourier transform, yielding a 
new method to compute a more accurate FFT. 
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1 Introduction 



A well-known fact is that the Fourier transform and some of its generalizations are fre- 
quently used tools in signal processing. The fractional Fourier transform PQ[2], which can 
be tracked down to the Wiener paper in 1929 [3] and is a particular case of the Linear 
Canonical Transform, which was derived in connection with canonical transformations 
in Quantum Mechanics [HE], has been object of renowned interest recently in the areas 
where the ordinary Fourier transform has been traditionally used. The work of Ozaktas 
et. al. [2] has established a standard framework to understand its properties and appli- 
cations. The techniques known collectively as time-frequency or space-frequency analysis 
take advantage of the intermediate representation of the signal and its Fourier transform 
to give a richer representation. In almost any domain where the usual Fourier transform 
is used, there is room for techniques based on the space-frequency analysis. As examples 
we have the solution of differential equations, quantum optics, signal processing, swept 
frequency and time varying filters, pattern recognition and the study of time-frequency 
distributions [5]-[TH]. Since this transform is a potentially useful tool for signal process- 
ing, the direct computation of the fractional Fourier transform in digital computers has 
become an important issue. In response to this need, an algorithm which computes a 
discrete fractional Fourier transform in 0(N log N) time has been published [6]. 
In this paper a discrete fractional Fourier transform is obtained in a vector space of di- 
mension iV by using some properties of the Hermite polynomials. The finite-dimensional 
vectors representing a Hermite function and its Fourier transform converge to their exact 
continuous values when N goes to infinity, therefore, a matrix operator can be proposed as 
a representation of the kernel of the fractional Fourier transform for functions other than 
the Hermite ones. This matrix, which plays the role of the kernel, yields a quadrature 
formula [HJ [H)J [TB] and a discrete form for the fractional Fourier transform. By using 
some asymptotic properties of the Hermite polynomials, this discrete fractional Fourier 
transform can be written in terms of the standard discrete Fourier transform through a 
diagonal congruence transformation. In this way, an efficient algorithm of 0(N log N) 
complexity for fast computations of fractional Fourier transforms is obtained. To dis- 
tinguish this discretization from other important contributions, particularly from those 
related with the discretization and development of fast algoritms to compute the linear 
canonical transformation [T7J [18], it is called XFT (extended Fourier Transform). The 
XFT can be evaluated in closed form in terms of exponentials for any complex value z of 
the unit circle \z\ < 1, though in the literature only values lying on the boundary \z\ = 1 
are usually considered. Due to the diagonal congruence transformations, the XFT can be 
understood complex- windowed FFT. 
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2 Quadrature of the Fractional Fourier Transform 



In this section, we reformulate the procedure followed in [TH[T5l[T6] to obtain a quadrature 
formula for the continuous fractional Fourier transform [I], yielding a discrete fractional 
Fourier transform, or as we call it, a extended fast Fourier Transform. 
Let us consider the family of Hermite polynomials H n (t), n = 0, 1, . . ., which satisfies the 
recurrence equation 

H n+1 (t) + 2n# n _i(*) = 2tH n (t), (1) 

with H-i(t) = 0. As it is well-known [19], from (pQ) follows the Christoffel-Darboux 
formula 
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Note that the recurrence equation ([T]) can be written as the eigenvalue problem 
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Let us now consider the eigenproblem associated to the principal submatrix of dimension 
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This is a general technique 
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to yield gaussian quadratures for orthogonal poly- 
nomials: the three-term recurrence equation is rewritten in matrix form to obtain or- 
thonormal vectors of M. N whose entries are given (in our case) in terms of the values of 
Hk(x), k = 0,1, . . . , N — 1, at the zeros of H^{x). We proceed by taking a similarity 
transformation to symmetrize Ti. Note that the diagonal matrix given by 



S = diag<J 1, -J=, . . . 
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generates the symmetric matrix 
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The recurrence equation ([I]) and formula ([2]) can be used to solve the eigenproblem 

Hu k = t k u k , k = 1,2, . . .,N, 

which is a finite-dimensional version of ([3]) . The eigenvalues t k are the zeros of (t) and 
the kth eigenvector u k is given by 



c k (s 1 H (t k ), s 2 Hi(t k ), s 3 H 2 (t h 



SNH N -i(t k )) J 



where Si, . . . , are the diagonal elements of S and c k is a normalization constant that 
can be determined from the condition u^u k = 1, i.e., from 

N-l 



2 sr^ H n (t k )H n (t k ) _ 
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Since H N {t k ) = and H' N {t k ) = 2NH N _i(t k ) ) the use of ([2]) yields 
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where we have used the fact that \H N _i(t k )\ = (—l) N+k H N _i(t k ). Thus, the components 
of the orthonormal vectors u k , k = 1, 2, . . . , N, are 



[u k 



-1) 



N+k 



! 2 N - n (N-l)\ H n ^(t h 



N(n-1)\ H N ^(t k) 



n 



1 N. 



(4) 



Let U be the orthogonal matrix whose kth column is u k and D(z) be the diagonal matrix 
D(z) = diagjl, z, z 2 , . . . , z 1 ^^ 1 }, where z is a complex number. Now, let us define the 
matrix 

T z = V2^U- l D(z)U 

whose components are given by 

N-l 



2 ) jk 



2tt 



(_l)J+fc 2 ^-i (JV-1)! 
N HN-i(tj)HN-i{th 
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5^ ^ni Hn ^ Hn ^- 



(5) 



This is the matrix representing the kernel of the fractional Fourier transform in a N- 
dimensional vector space, as we show next. 
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Asymptotic formulae 



Let us look for the asymptotic form of the components (jSJ) of T z . First note that the 
asymptotic expression for -£/jv(t) in the oscillatory region is (Eq. (8.22.8) of [19|) 

H N (t) = r(N/2 + l) e * /2 {™<V2N + 1 t - iVvr/2) + 0(N^ 2 )) . (6) 

This gives an approximate form for the zeros of H^(t), i.e., 

(2k- N - 1\ 7T 
- 77, (7) 



f 2N 

k = 1, 2, . . . , N. Thus, the use of © and © yields 
Therefore, for N great enough, ([5]) can be written as 

lAT-l rr/JV+l\n2 °° 



^ ' n=0 

and finally, Stirling's formula and Mehler's formula [22] produce the result 



" V GXP I 2(1^2) J fe ' (8) 

where Ai^ is the difference between two consecutive asymptotic Hermite zeros, i.e., 

At k = t k+ i —t k = _ - . (9) 

Let us consider now a complex- valued function defined for t 6 M and let us form the 
vector 

9 = {g{ti),g{t 2 ),...,g{t N )) T . 

Therefore, the multiplication of the matrix T z by the vector g gives the vector G with 
entries 

r VV^ r^V- f (l + ^)(^ + ^)-4ti^ ..... nm 
g j = zJFz)jk9(tk) ^ y ^—^ Z^ ex P ( 2(1 - z 2 ) j 9 ^ Atk > ( 10 ) 

k=i v fc=i V ^ ' / 

for j = 1, 2 . . . , N. Note that this equation is a Riemann sum for the integral 

^r,n , /"°° ( (1 + z 2 )(t 2 + t' 2 ) -Att' z\ /A . , , x 

W^ =Ji ?/ ex P 1 )9{t')dt', \z\<\, (11) 
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that is, 

N 

& x \gtf)M - ^2(Fz)jkg(t k ), N^oo. (12) 

k=l 

Note that JP z [g(t'),t] is the continuous fractional Fourier transform pQ of g(t') up to 
a constant. Thus, the matrix T z is a discrete fractional Fourier transform. Since the 
fractional transform is an extension for the Fourier transform, T z is termed XFT. It 
should be noted that the XFT is a discretization of the fractional Fourier transform for 
any complex value z of the unitary circle | z | < 1 and not only for z lying on its boundary, 
as it is usually considered. Note that for the XFT, the argument <p of z = rexp(i(p) is 
real and not complex as it is used in some applications |23j . 
In the case z = ±i, ([8]) becomes 

F jk = (F ±i ) jk ~ e ±lt ^At k , 

that is, the XFT becomes a discrete Fourier transform. It has been previously obtained 
in [13] where some numerical examples are given, and it has been applied to the analysis 
of brain signals [24J . 



3 The fast XFT 

The FFT can be used to obtain a fast algorithm for the XFT as following. Since the matrix 
T z represents a quadrature for the fractional Fourier transform expected to converge when 
a great number of nodes are used, we consider the asymptotic form ([S]) which can be 
written in more detail as 



{Fz)jk = \l 1 ~ z2 exp(-/itj) exp(utjt k ) exp(-/it fc )Ai fe , (13) 
where we have used the definitions 

1 + z2 2z 
2(1 — z z ) 1 — z z 

Note that we have replaced the approximately equal sign "~" by the equal sign "=" in 
f fT3|) redefining T z . In order to show the main differences of the XFT with the usual FFT, 
we consider first the case of the standard Fourier transform. 



3.1 XFT as an improvement of the FFT 



As noted above, the case z = i in f|T3|) corresponds to a discrete Fourier transform 
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(15) 



where now j, k — 0, 1, 2, . . . , N — 1, and we have used ([7]) and Q). Since X^feLi Fjk9{tk) is 
a quadrature and therefore, an approximation [cf. Eq. [T2] of 



we can use the basic property of the Fourier transform 
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to rewrite (11 5p in the more convenient form 
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(16) 



(17) 



in the understanding that this matrix yields a scaled Fourier transform. If we choose 
a = A/tt, Eq. (iTTj) can be written as 



7re 2 



4/ttM 



where j, k — 0, 1, 2, . . . , N — 1. In matrix form, and defining F = i 7 ^/^ to simplify the 
notation, 



ire 2 at 
\/2iV 



-5 , J D F S', 



18) 



where 5 is the diagonal matrix whose nonzero elements are exp(— nr(iV — l)j/N), j = 
0, 1, . . . , N — 1 and D F (-) is just the FFT of the argument. 

Note that a very important contribution of our result is that it truly enables the extension 
of the FFT algorithms to a more general case. The fact that computing F(-) is essentially 
computing the FFT of the vector S(-), means that almost at no cost, modern instruments 
would be able to increase ther capabilities for digital processing. 
The inverse of F is easily obtained: 



7T 



2tt 



N-l 



k 



N 



where j,k = 0, 1, 2, . . . , iV - 1. 

In the applications, we have to remind that F gives a scaled transform. The following 
simple algorithm incorporate these ideas. 
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Algorithm 1 



To compute an approximation G = (G\, G2, ■ ■ ■ , Gn) t of the Fourier transform 
ffT6l) of the vector 

9 = (9i,92, ■ ■ -,9n) T - 

1. For j, fc = 0,l,2,...,iV-l, 

(a) Compute (D F )j k = e l ^wi k (the discrete Fourier transform). 

(b) Compute the diagonal matrix S according to Sjk = e~ mE N^^5jf.- 

2. Obtain the approximation Gj to G(-tj) by computing the matrix-vector 
product 

(JV-1)2 

ire 2 n 

G = -j= — SDp(Sg), (19) 

V2N 

with a standard FFT algorithm. 



Either if the input vector g is given by the values of a function g(t) at t k , or not, we can 
represent the approximation to the Fourier transform G(u), as given in a plot of points 



In the examples given below we plot the real and imaginary parts of the vector G and the 
exact transform G(-uS). We give first two examples of non-periodic/singular functions 
|25j . To compare the performance of XFT with that of FFT, we show the plot of the 
corresponding FFT in the first example. 



Example 1 

Consider the pair of Fourier transforms 

g(t) = cos(t 2 ), G(uo) = v^rcosl 



0J 2 — 7T, 



4 



Figures 1 and 2 shows the numerical convergence attained by the XFT for N=512 and 
1024 respectively and it is compared with the FFT output. 
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Figure 1: A Part of the plots corresponding to (tj, G(-tj)) (solid line) and (tj, Gj) (dashed 
line). The max-norm of the overall error is 2.11 for N = 512. B The standard FFT 
computed for the same function and number of points. 
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Figure 2: A Part of the plots corresponding to (tj, G(-tj)) (solid line) and (tj, Gj) (dashed 
line). The max-norm of the overall error is 2.08 for N = 1024. B The standard FFT 
computed for the same function and number of points. 



Example 2 



As an example of a singular function, let us consider the pair of Fourier transforms (the 
integral is a Cauchy Principal Value) 

Figures 3 and 4 shows the good performance of the XFT for N=512 and 1024 respectively. 
In this cases we have small absolute errors. 
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Figure 3: A Real part and B imaginary part of the continuous Fourier transform (solid 
line) compared with the XFT (dashed line) computed with N = 512. The function f(t) 
is that given in Example 2. The max-norm of the overall error is 0.4262 for both the real 
and imaginary part. 





Figure 4: A Real part and B imaginary part of the continuous Fourier transform (solid 
line) compared with the XFT (dashed line) computed with N = 1024. The function f(t) 
is that given in Example 2. The max-norm of the overall error is 0.105 for both the real 
and imaginary part. 



Example 3 

This example is concerned with the performance of the XFT acting on sinusoidal wave- 
forms of integer-harmonic frequency . Our starting point is again Eq. (|17p . By changing 
the order of the indexes, (IP7|) can be written as 

{F)it = 7m e 



where 

'0,±1,±2, ••• ,±(JV-l)/2, Nodd, 
±1/2, ±3/2,- •• ,±{N- l)/2, Neven. 
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Let N be odd, and assume that the signal to be transformed is an integer-harmonic signal. 
For instance, 

g k = cos(fcr m ), r m = ——, (20) 
where k = 0, ±1, ±2, . . . , ±(N - 1) /2. Then, the formulae [26] 

, sm(Nx/2)cos((N + l)x/2) sin(iVx/2) sin((iV + l)x/2) 

> cos fei = ; — — , > sin b = ; — — , 

t! y ' sin(x/2) sin(x/2) 

can be used to prove that the XFT transform of (120]) is 

{N - 1)/2 „ [W 

Gj= ^2 (F) jk cos(kT m ) = Ti J —{5j- m + 5 jm ), (21) 

fe=-(iV-l)/2 

where j = 0, ±1, ±2, ■ ■ • , ±(iV— 1)/2, corresponding to two pulses centered at the integers 
m and — m respectively. 

Now assume that N is even and take a half- integer harmonic signal, i.e., one of the same 
form as ( 1201) but with k = ±1/2, ±3/2, . . . , ±(iV — l)/2. A similar calculation as above 
shows that Gj has the same form of ( 12TI) but the sum runs over half-integers and the 
pulses are now centered at the half-integers m and — m respectively. 
In the two above cases the dispersion is zero, but this does not happen if we take iV odd 
and consider a half-integer harmonic signal, or iV even and consider an integer harmonic 
signal. Finally note that these processes can be reversed: the XFT of the corresponding 
pulses yields integer or half-integer harmonic signals. 

We give another numerical case with a cosine waveform with an arbitrary frequency. Let 
us consider the example g(t) = cos(5.156t). In Fig. 5 we show the XFT output for (a) 
iV = 1024 and (b) N = 2048. As it can be seen, the leakage practically vanishes for the 
latter case. As a measure of the leakage we give the mean value 



1 N 



k=l 



where the indexes m and m! correspond to the frequencies where |G| attains their maxima. 
Thus, we find that /i = 0.14105 for N = 1024 and /i = 0.00276 for N = 1024. 
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Figure 5: Fourier transform of cos(5.156£). The maximum occurs at u — 5.17072 for A 
N = 1024 and at tu = 5.15625 for B iV = 2048. 



3.2 XFT as a fast discrete fractional Fourier transform 

To obtain a discrete and fast implementation of the continuous fractional Fourier trans- 
form, we follow the same procedure as above. The definition of the fractional transform 
given by ffTTT) is \/2ti times the one given in [TJ. This has to be taken into account for 
comparing purposes. 

Our starting point is the definition of the fractional transform given in ( TTTT) 



/oo 
e«* 'e~^ 'g(t')dt> \ 

where we have used the definitions (fl4"l) . The simple scaling 



G *( at ) = I e avtt 'e-^' 2 g(t')dt', (22) 

and the discretization of the kernel 



(Fz)jk = \/ ?~~2 exp(- y ua 2 t|)exp(az/t i t fc )exp(- y u^)At A ., 

allows us to implement a fast algorithm to compute the XFT in terms of the FFT. This 
can be done by choosing a = i2(l — z 2 )/(ttz), because then 



2 2i 

(Fz)jk = \l ; o ex P(-/ ia2t ?)( F )ifc ex P(-^D ? a = — (! - ^ 2 )> 

V 1 — z z J nz 

where F, given by fflBl . is the fast XFT for the case of the standard Fourier transform. 
In order to simplify the notation, let T z be j^ 2l ^ z )/( 7r2 )_ Since ^2^ =1 (^ z )jkf(tk) is a 
quadrature of fl22|) with a = i2(l — z 2 )/(tcz) [cf. Eq. [12], T z is a fast discrete fractional 
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Fourier transform which gives an approximation for the scaled function G z (at) at the 
nodes t k . Note that the inverse transform can be obtained in terms of F~ l . These ideas 
are incorporated in the following 



Algorithm 2 



To compute an approximation G z = (G z ±, G Z 2, ■ ■ ■ , G z n) t of the fractional Fourier 
transform JP z [g(t'), i2(l — z 2 ) / (j[z)tj\ of the vector 

9 = (9i,92, ■ ■ -,9n) T , 
at the complex value z, \z\ < 1. 

1. For given N and z, set 

(a) /i = (1 + z 2 )/[2(l - z% a = 2i(l - z 2 )/(ttz). 

(b) t k = n(2k- N -1)/[2(V2N)}, k = 1, 2, . . . , JV 

2. For j,fc = 0,l,2,...,iV-l, 

(a) Compute (Dp)jk = e l ^^ k (the discrete Fourier transform). 

(b) Compute the diagonal matrices S\ and S2 according to 

(S,)^ = e-^~ l ^5 jk , (S 2 ) jk = e-^^S jk . 



3. Obtain the approximation G Z j to G z (-tj) by computing the matrix-vector 
product 



2 tip' 1 2 N 

O, - xIy—; -= SMSig), (23) 



with a standard FFT algorithm. 



We first test this fast transform on two well-known examples for which \z\ = 1 (see Ref. 

DO)- 



Example 4 

Consider the pair of fractional Fourier transforms 

g{t) = exp(-t 2 /2 + pt), G z (t) = exp(-t 2 /2 - (i/2)/3 2 e itp sin <p + /3te i(fi ), z = e itfi '. 

Fig. 6 shows a part of the plot of the real and imaginary parts of the XFT compared with 
G z (atj) for this case. The max-norm of the error is of order 10~ 12 for both cases. 
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Figure 6: A Real part and B imaginary part of the fractional Fourier transform (solid 
line) compared with the XFT (dashed line) computed with N = 512. The function f(t) 
is that given in Example 4 for (3 = 2. 



Example 5 

Consider the pair of fractional Fourier transforms 

g(t) = 1, G z {t) = exp[i(t 2 tan^ - ip)/2]/y/cos(p, z = e i(p . 

Fig. 7 shows the whole plot of the real and imaginary parts of the XFT compared with 
G z (atj) for this case. The max-norm of the overall error is 1.3282 for the real part and 
1.42694 for the imaginary part. 




Figure 7: A Real part and B imaginary part of the fractional Fourier transform (solid 
line) compared with the XFT (dashed line) computed with N = 512. The function f(t) 
is that given in Example 5. 
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Example 6 



This example consider the fractional Fourier transform of the rectangle function. The 
cases correspond to the fractional Fourier transform as defined in the literature PQ, for 
the values of the argument (p given in the caption. This figure shows how the Fourier 
transform evolves as the argument of z — exp(i(p) takes decreasing values. 
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Figure 8: Fractional Fourier transform as given in ffTTT) for the rectangle function rect(t), 
computed with N = 512. We plot in A the exact Fourier transform (solid line) against 
the XFT (dashed line). In this case we have <p = tt/2. In the cases B, C and D, the 
argument ip is 1,1/2 and 1/4 respectively. Note that the spectrum peaks becomes greater 
as the argument (p decreases. 



4 Conclusions 

As it is known, there are not universal methods. Some points for and against the fast 
algorithm to compute the fractional Fourier transform presented in this paper have to be 
noticed. First, note that the XFT and their inverse have simple forms given by explicit 
matrix elements that can be computed easily for complex values of the parameter z = 
rexp(i<p) lying inside the unitary circle, not only on the boundary. On the other hand, the 
XFT can be interpreted as a complex-windowed FFT, as fast as the latter and certainly, 
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it is an improvement of the FFT since it performs pretty well even with singular functions 
or non-periodic functions and detects the well-known integer-harmonic spectrum as the 
FFT does. Besides, it is able to detect half-integer harmonic frequencies with no leakage. 
However, the convergence of the XFT to the fractional Fourier transform can be assured 
only if the function to be transformed satisfies the quadrature formula (fl2l) . Besides, 
the XFT also presents some pitfalls and artifacts as the FFT: it presents leakage for 
some sinusoidal signals and yields spurious frequencies if the input signal has very high 
frequencies and the number of samples is small. 
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